Resumen

La accidentalidad vial es un fenómeno que se está convirtiendo en un problema de salud pública dado que toca directamente al ser humano y detrás de cada cifra hay un drama familiar, por esta razón las entidades gubernamentales tienen que ser más ambiciosos en materia de campañas, entendiendo que hay una gran inversión desde el orden nacional para resolver esta problemática, pero no es suficiente, por lo cual se debe insistir más en el tema de seguridad vial desde lo local. Al respecto, la alcaldía de Medellín a venido recolectando datos sobre accidentalidad vehicular en el municipio desde el año 2014 con el fin de describir el comportamiento y los lugares más afectados por este problema. En este trabajo se utiliza información de accidentes de 2014 a 2018 para describir estadisticamente los accidentes presentados y tratar de predecir accidentes presentados en las diferentes epocas del año y sitios del municipio.

  1. Introducción

La Estadística se ocupa del manejo de la información que pueda ser cuantificada. Implica esto la descripción de conjuntos de datos y la inferencia a partir de la información recolectada de un fenómeno de interés. La función principal de la estadística abarca: Resumir, Simplificar, Comparar, Relacionar, Proyectar. Constantemente se generan datos que se alojan en los sistemas de información y bases de datos públicas y privadas. La Alcaldía de Medellín cuenta con una página web donde alojan bases de datos públicas en diferentes temas sobre la ciudad. El objetivo de este trabajo es utilizar la Estadística y el desarrollo computacional para analizar los datos de incidentes georreferenciados de 2014 a 2018 de la ciudad de Medellín para predecir accidentes del año 2019 con el fin de descubrir factores que se relacionan con los accidentes y ayudar de alguna manera a las entidades públicas a tomar decisiones en temas de seguridad vial.

Documentación Soporte

  • En el enlace se encuentra alojado el aplicativo web VehApp que permite visualizar la accidentalidad y pronosticar Enlace: Aplicacion web VehApp

  • En el enlace se encuentra alojado el Github principal que contiene todos los archivos y scripts utilizados para la construcción e implementación de este trabajo Enlace: Github principal

  1. Metodología

En inteligencia artificial un problema se define como la existencia de discrepancias entre un estado actual y un estado objetivo. La resolución recurrente y especializada de ciertos tipos de problemas permite identificarlos en forma de tareas. Para abordar este problema se usa la metodología para la Analítica de datos (ASUM-DM) ya que establece una forma probada para llevar un problema a una solución de analítica. ASUM-DM de IBM se basa en la metodología CRISP-DM y hace enfasis en las nuevas prácticas en la ciencia de datos como el uso de volumenes de datos muy grandes, la incorporación de análisis de texto, en el modelado predictivo y la automatización de algunos procesos. El ciclo de elaboración se basa en 9 principales pasos:

  1. Entendimiento del negocio
  2. Descubrimiento y Entendimiento de datos
  3. Preparación y gestión de datos
  4. Validación de datos
  5. Análisis
  6. Prototipo
  7. Validación y evaluación
  8. Nueva adquisición de datos
  9. Ciclo de publicación

  1. Desarrollo

3.1 Entendimiento del negocio

Más de un millón de personas muere cada año en accidentes automovilísticos, Estos incidentes ocurren por una cantidad de razones, pero el error humano está en el centro. Los conductores manejan más rápido de lo que deberían en ciertas condiciones atmosféricas, toman decisiones equivocadas y no logran anticipar riesgos potenciales. El conocimiento de las variables involucradas en los accidentes de tránsito pueden facilitar la formulación de políticas específicas para reducir la creciente accidentalidad vial. Asi, el objetivo de este trabajo es describir las características de tiempo, lugar, vehiculo y demás variables que puedan inferir en los accidentes de tránsito de la ciudad de Medellín entre 2014 y 2018, buscando aportar elementos al conocimiento de este problema que sean útiles para quienes allí formulan las políticas de prevención y control de morbilidad y mortalidad.

3.2 Descubrimiento y entendimiento de datos

La Alcaldía de Medellín en el portal GeoMedellín suministra las bases de datos anuales para incidentes de tránsito registrados por la Secretaria de Movilidad. Los datos se encuentran en este enlace.. Para el entendimiento del conjunto de datos se define el accidente de transito como: “evento, generalmente involuntario, generado al menos por un vehiculo en movimiento, que causa daños a personas y bienes involucrados en él e igualmente afecta la normal circulación de los vehiculos que se movilizan por la via o vias comprendidas en el lugar o dentro de la zona de influencia del hecho”

La estructura de las variables que contiene los datos es la siguiente:

Variable Descripción Tipo Detalle
X Coordenada X (longitud) de la ubicación del accidente float -
Y Coordenada Y (latitud) de la ubicación del accidente) float -
OBJECTID Identificación del registro integer -
RADICADO Identificación única del accidente ante la Secretaría de Movilidad string -
FECHA Fecha de la ocurrencia del accidente string -
HORA Hora aproximada del accidente datetime -
DIA Número día dentro del mes del accidente integer -
PERIODO Año de la ocurrencia del accidente integer -
CLASE Tipo de accidente String Clases: atropello, caída del ocupante, choque, incendio, volcamiento, otro
DIRECCION Dirección detallada de la ubicación del accidente String -
DIRECCION_ENC Dirección encasillada de la ubicación del accidente string Formato único de direcciones de la Alcaldía de Medellín
CBML Identificación única del lote más cercano a la ubicación del accidente string Acrónimo de comuna, barrio, manzana, lote
TIPO_GEOCOD Tipo de ubicación catastral string Formato de la oficina de catastro de la Alcaldía de Medellín
GRAVEDAD Gravedad del accidente string Clases: herido, muerto, solo daños
BARRIO Barrio donde ocurrió el accidente string 340 niveles dentro de la variable
COMUNA Comuna donde ocurrió el accidente string 85 niveles dentro de la variable
DISENO Tipo de diseño en donde ocurrió el accidente string Clases: ciclo ruta, glorieta, intersección, lote o predio, paso a nivel, paso elevado, paso inferior, pontón, puente, tramo de vía, túnel, vía peatonal
DIA_NOMBRE Nombre del día de la semana del accidente string -
MES Número de mes dentro del año del accidente integer -
MES_NOMBRE Nombre del mes del accidente string -
X_MAGNAMED - - -
Y_MAGNAMED - - -
LONGITUD Coordenada longitud de la ubicación del accidente float -
LATITUD Coordenada latitud de la ubicación del accidente float -

3.3 Preparación y gestión de datos

El conjunto de datos llamado medellin cuenta con 228.693 registros y 24 variables que corresponden a accidentes desde el 2014 a 2018. Para el inicio de la preparación de ellos se analizará que variables no hay necesidad de tener en cuenta, cuantas variables hay que transformar y que otras se pueden agregar. Para esto se hace uso de Tidyverse en el lenguaje R para la limpieza de datos y otros paquetes necesarios.

Desde el portal de GeoMedellín se descargan datos de accidentalidad por año desde 2014 a 2018 en archivos csv, se ingresan a una carpeta del computador y se hace lectura de cada uno de ellos luego se juntan para crear solo un conjunto de datos que abarca los 5 años.

# Cargar libreria
library(data.table);library(purrr);library(dplyr)        
library(plotly);library(tidyr);library(stringr)      
library(lubridate);library(sf);library(knitr);library(kableExtra)

## lectura de datos

# lista archivos
setwd("Datos")
archivos_datos <- list.files(pattern = ".csv", include.dirs = T, recursive = F)

# lectura de los archivos (2014 a 2018)
datos <- map(archivos_datos, fread, sep = ",", encoding = "UTF-8", colClasses = "c")

# juntar por filas para crear un solo conjunto de datos
medellin <- bind_rows(datos)

# ver datos
kable(head(medellin), align = "ccccc") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "300px")
X Y OBJECTID RADICADO FECHA HORA DIA PERIODO CLASE DIRECCION DIRECCION_ENC CBML TIPO_GEOCOD GRAVEDAD BARRIO COMUNA DISENO DIA_NOMBRE MES MES_NOMBRE X_MAGNAMED Y_MAGNAMED LONGITUD LATITUD
-75.6027257778941 6.21901553077807 211279 1423828 2014/01/01 00:00:00+00 02:20:00 01 2014 Choque CR 80 CL 8 CR 080 008 000 00000 1611 Malla vial aproximada: CR 81-7 HERIDO Loma de los Bernal Belén Tramo de via MIÉRCOLES 01 831190.08 1179689.5 -75.60272578 6.21901553
-75.5681847049393 6.26000899508826 211280 1423839 2014/01/01 00:00:00+00 12:50:00 01 2014 Atropello CR 53 CL 61 CR 053 061 000 00000 1003 Malla vial HERIDO Jesús Nazareno La Candelaria Tramo de via MIÉRCOLES 01 835013.63 1184223.79 -75.56818471 6.260009
-75.5499431117657 6.26476471445961 211281 1423840 2014/01/01 00:00:00+00 01:00:00 01 2014 Atropello CR 39 CL 70 CR 039 070 000 00000 0308 Malla vial HERIDO Manrique Oriental Manrique Tramo de via MIÉRCOLES 01 837032.71 1184749.86 -75.54994311 6.26476471
-75.6076107982279 6.23432736551888 211282 1423849 2014/01/01 00:00:00+00 12:37:00 01 2014 Atropello CL 32 CR 84 CL 032 084 000 00000 1617 Malla vial HERIDO Las Mercedes Belén Tramo de via MIÉRCOLES 01 830649.47 1181383.24 -75.6076108 6.23432737
-75.5796854467144 6.29996777075499 211283 1423890 2014/01/01 00:00:00+00 10:40:00 01 2014 Caída de Ocupante CR 80 CL 101 CR 080 101 000 00000 0603 Malla vial HERIDO Doce de Octubre No.2 Doce de Octubre Tramo de via MIÉRCOLES 01 833740.78 1188643.8 -75.57968545 6.29996777
-75.5554293007481 6.28449892412236 211284 1423892 2014/01/01 00:00:00+00 04:00:00 01 2014 Choque CR 48 CL 93 CR 048 093 000 00000 0401 Malla vial HERIDO Berlin Aranjuez Interseccion MIÉRCOLES 01 836425.43 1186932.72 -75.5554293 6.28449892

Con la descripción de las variables que contienen los datos se identifica las siguientes correcciones que se deben hacer a los datos:

  • Convertir la FECHA en formato YYYY-MM-DD

  • Eliminar datos vacios en la variable CLASE, quitar tildes y unificar nombres de niveles

  • Los datos vacios en la variable DISENO se cambian por OTRO, quitar tildes y unificar nombres de niveles

  • Quitar tildes a nombres de dias

# Organizar fecha
medellin$FECHA <- as.Date(medellin$FECHA) # Formato YYYY-MM-DD

# Eliminar datos vacios de la variable CLASE y correccion de niveles
medellin <- medellin[-which(medellin$CLASE == ""),]

# Correccion tildes
medellin$CLASE <- iconv(medellin$CLASE, from = "UTF-8", to = "ASCII//TRANSLIT")

# Quitar "DE" (nivel "Caída de ocupante" variable CLASE)
medellin$CLASE <- gsub(pattern = "Caida de Ocupante", replacement = "Caida Ocupante", x = medellin$CLASE) 

#Cambiar nivel vacio por otro en variable DISENO
medellin$DISENO <- ifelse(medellin$DISENO == "","Otro",medellin$DISENO) 

# Quitar tildes en niveles de la variable DISENO
medellin$DISENO <- iconv(medellin$DISENO, from = "UTF-8",to="ASCII//TRANSLIT") #Quitar tildes

# Quitar tildes a nombres de dias en variable DIA
medellin$DIA_NOMBRE <- iconv(medellin$DIA_NOMBRE, from = "UTF-8", to = "ASCII//TRANSLIT")

# ver estructura variables
str(medellin)
## Classes 'data.table' and 'data.frame':   228687 obs. of  24 variables:
##  $ X            : chr  "-75.6027257778941" "-75.5681847049393" "-75.5499431117657" "-75.6076107982279" ...
##  $ Y            : chr  "6.21901553077807" "6.26000899508826" "6.26476471445961" "6.23432736551888" ...
##  $ OBJECTID     : chr  "211279" "211280" "211281" "211282" ...
##  $ RADICADO     : chr  "1423828" "1423839" "1423840" "1423849" ...
##  $ FECHA        : Date, format: "2014-01-01" "2014-01-01" ...
##  $ HORA         : chr  "02:20:00" "12:50:00" "01:00:00" "12:37:00" ...
##  $ DIA          : chr  "01" "01" "01" "01" ...
##  $ PERIODO      : chr  "2014" "2014" "2014" "2014" ...
##  $ CLASE        : chr  "Choque" "Atropello" "Atropello" "Atropello" ...
##  $ DIRECCION    : chr  "CR 80 CL 8" "CR 53 CL 61" "CR 39 CL 70" "CL 32 CR 84" ...
##  $ DIRECCION_ENC: chr  "CR  080   008  000 00000" "CR  053   061  000 00000" "CR  039   070  000 00000" "CL  032   084  000 00000" ...
##  $ CBML         : chr  "1611" "1003" "0308" "1617" ...
##  $ TIPO_GEOCOD  : chr  "Malla vial aproximada: CR 81-7" "Malla vial" "Malla vial" "Malla vial" ...
##  $ GRAVEDAD     : chr  "HERIDO" "HERIDO" "HERIDO" "HERIDO" ...
##  $ BARRIO       : chr  "Loma de los Bernal" "Jesús Nazareno" "Manrique Oriental" "Las Mercedes" ...
##  $ COMUNA       : chr  "Belén" "La Candelaria" "Manrique" "Belén" ...
##  $ DISENO       : chr  "Tramo de via" "Tramo de via" "Tramo de via" "Tramo de via" ...
##  $ DIA_NOMBRE   : chr  "MIERCOLES" "MIERCOLES" "MIERCOLES" "MIERCOLES" ...
##  $ MES          : chr  "01" "01" "01" "01" ...
##  $ MES_NOMBRE   : chr  "" "" "" "" ...
##  $ X_MAGNAMED   : chr  "831190.08" "835013.63" "837032.71" "830649.47" ...
##  $ Y_MAGNAMED   : chr  "1179689.5" "1184223.79" "1184749.86" "1181383.24" ...
##  $ LONGITUD     : chr  "-75.60272578" "-75.56818471" "-75.54994311" "-75.6076108" ...
##  $ LATITUD      : chr  "6.21901553" "6.260009" "6.26476471" "6.23432737" ...
##  - attr(*, ".internal.selfref")=<externalptr>
Variables BARRIO y COMUNA

La variable barrio y comuna es muy importante para el objetivo, por eso se hace un análisis completo de como están estructurados sus niveles y las reformas que se le deben hacer. Se encontró que la variable COMUNA tiene nombres erroneos y además contiene 19.736 registros vacios. En la variable BARRIO se tienen 19.766 registros vacios y algunos de sus niveles contienen numeros que no indican bien el nombre de la ubicación. Para la solución de este problema se usó la variable CBML que es una identificación única del lote más cercano al accidente conformado por un Acrónimo de comuna, barrio, manzana, lote, con esta variable se puede extraer un registro con los 4 primeros números para obtener la comuna y el barrio. La oficina de catastro cuenta con un shapefile donde se registran verdaderamente los nombres de barrios y comunas identificados con su codigo que en los datos de accidentes serían esos 4 primeros digitos del CBML, así que se hace un cruce entre los datos de accidentes y los datos de catastro uniendo por el código. Como observación la variable CBML cuenta con 19.757 datos vacios y 185 registros en los cuales no cuenta con un registro CBML claro, que ayude a identificar el barrio y la comuna.

Después de tener organizados barrios y comunas, se eliminan las variables que no son necesarias como LATITUD Y LONGITUD dado que X y Y representan lo mismo, RADICADO, CBML, TIPO GEOCOD, X_MAGNAMED, Y_MAGNAMED, DIRECCION, DIRECCION_ENC y MES_NOMBRE también se eliminan y se realiza un ajuste en el nombre de las variables. El nuevo conjunto de datos queda compuesto de 208.745 registros y 15 variables

#==========================
# ARREGLOS COMUNA Y BARRIO
#===========================
setwd("Datos")
catastro <- read_sf("Limite_Barrio_Vereda_Catastral.shp")

# Crear un codigo desde CBML que deje solo comuna y barrio
medellin <- mutate(medellin, CB = str_sub(CBML,1,4))

# ingresar las varaibles de identificacion de lugar dadas por catastro en el archivo de datos
medellin <- inner_join(medellin, select(catastro, CODIGO, NOMBRE_COM, NOMBRE_BAR),
                        by = c("CB" = "CODIGO"))

# Quitar duplicados por el inner_join. Pasamos de 228687 registros a 208745 registros
medellin <- medellin[!duplicated(medellin[,c("OBJECTID")])]

# Quitar tildes y volver a mayuscula los nombres de barrios
medellin$NOMBRE_BAR <- iconv(toupper(medellin$NOMBRE_BAR), from = "UTF-8", to = "ASCII//TRANSLIT")

# Elegimos las columnas que necesarias
medellin <- medellin[,c("X","Y","OBJECTID","FECHA","HORA","DIA","PERIODO","CLASE","GRAVEDAD",
                        "DISENO","DIA_NOMBRE","MES","CB","NOMBRE_COM","NOMBRE_BAR")]

# Cambiemos algunos nombres de las variables
names(medellin)[c(1,2,14,15)] <- c("LONGITUD","LATITUD","COMUNA","BARRIO")

# Asi quedan los datos
head(medellin)
Nuevas variables

En las nuevas variables a ingresar se considera días especiales que corresponden a fechas que se caracterizan por ser grandes eventos en los que la comunidad participa de manera completa o mayoritaria y que de alguna manera afecta la movilidad, para esto se crea un archivo plano csv que tiene el registro de cuales han sido fechas especiales desde 2014 a 2018 y el tipo de festividad, para esta variable se tuvo en cuenta las siguientes fechas especiales.

  • BRUJOS: 31-oct día Halloween
  • FESTIVO: Días festivos colombianos
  • MADRE: Día de la Madre
  • NAVIDAD: 24 y 31 de Diciembre
  • SANTA: Todos los días que conforman la Semana Santa
  • VACACIONES: Semana de Octubre donde estudiantes de colegios salen a vacaciones

Otra variable nueva a ingresar se crea a partir de la hora del accidente. En la ciudad de Medellín se considera como hora a pico al tiempo en que la mayoría de los habitantes se dirigen al trabajo, estudio o de regreso a casa, asi el intervalo de tiempo que se toma como hora pico es de 5am a 8am y de 5pm a 8pm. Cabe resaltar en la estructuración de la variable hora no está especificado a que parte del día ocurrió el accidente, es decir, con la hora no se puede saber si ocurrió en la mañana, en la tarde o en la noche.

Otras dos variables nuevas a ingresar corresponden al estado del clima:

  • HUMEDAD_REL: Humedad relativa - Es una relación, expresada en porcentaje, de la cantidad de humedad atmosférica presente relativa a la cantidad que estaría presente si el aire estuviera saturado, esta relación es determinante para que se produzca evaporación o precipitación (lluvia)
  • PLUVIO: Pluviometria [mm] - Mide la cantidad de precipitaciones (lluvias) caídas en un lugar durante un tiempo determinado

Se debe advertir que la variable HUMEDAD_REL tiene registros desde 2014-12-04 a 2018-12-01 y la variable PLUVIO desde 2014-01-01 a 2018-12-01

Los nuevos datos quedan conformados por 208.745 observaciones y 19 variables. Para finalizar con el arreglo de los datos se realiza una estructuracion en el tipo de variable necesaria para los análisis.

#==============================
# Ingreso de fechas especiales
#==============================
setwd("Datos")
# Lectura registros de fechas especiales desde el 2014 a 2019
especiales <- read.csv("Otros/Fechas_especiales.csv", sep = ";", header = T)

# Convertir en formato YYYY-MM -DD
especiales$FECHA <- as.Date(especiales$FECHA, format = "%d/%m/%Y")

# unir las fechas especiales a medellin, los que nos son fechas especiales los pone como NA
medellin <- merge(x = medellin, y = especiales, by = "FECHA", all.x = T)
medellin$ESPECIAL <- ifelse(is.na(medellin$ESPECIAL),"NO","SI")

# Variable HORA_PICO
medellin$HORA <- hms(medellin$HORA)
medellin$HORA_PICO <- ifelse(hour(medellin$HORA) %in% c(5,6,7), "SI","NO")


# Variable HUMEDAD_REL y PLUVO

# Humedad Relativa
HUMEDAD <- read.csv('Otros/Datos IDEAM/HumedadRelativa.csv.csv')
HUMEDAD <- data.frame(FECHA=as.Date(HUMEDAD$Fecha),HUMEDAD_REL=HUMEDAD$Valor)

# Pluviometria
PLUV <- read.csv('Otros/Datos IDEAM/Pluviometria.csv.csv')
PLUV <- data.frame(FECHA=as.Date(PLUV$Fecha),PLUVIO = PLUV$Valor)
PLUV <- PLUV[!duplicated(PLUV[,c("FECHA")]),]

medellin <- merge(x = medellin, y = HUMEDAD, by = "FECHA", all.x = TRUE)
medellin <- merge(x = medellin, y = PLUV, by = "FECHA", all.x = TRUE)

# Tipo variables
medellin$LONGITUD <- as.numeric(medellin$LONGITUD)
medellin$LATITUD <- as.numeric(medellin$LATITUD)
medellin$DIA <- as.integer(medellin$DIA)
medellin$PERIODO <- as.integer(medellin$PERIODO)
medellin$CLASE <- as.factor(medellin$CLASE)
medellin$GRAVEDAD <- as.factor(medellin$GRAVEDAD)
medellin$DISENO <- as.factor(medellin$DISENO)
medellin$DIA_NOMBRE <- as.factor(medellin$DIA_NOMBRE)
medellin$MES <- as.integer(medellin$MES)
medellin$COMUNA <- as.factor(medellin$COMUNA)
medellin$BARRIO <- as.factor(medellin$BARRIO)
medellin$ESPECIAL <- as.factor(medellin$ESPECIAL)
medellin$HORA_PICO <- as.factor(medellin$HORA_PICO)

# Visualizar datos
head(medellin)

Y su nueva estructura de variables es la siguiente:

Variable Descripción Tipo Detalle
FECHA Fecha ocurrencia accidente Date -
LONGITUD Coordenada longitud ubicación accidente integer -
LATITUD Coordenada latitud ubicación accidente integer -
OBJECTID Identificación única del registro string -
HORA Hora aproximada del accidente datetime -
PERIODO Año de la ocurrencia del accidente integer -
CLASE Tipo de accidente FACTOR 7 niveles: Atropello, Caida Ocupante, Choque, Choque y Atropello, Incendio, Volcamiento, Otro
GRAVEDAD Tipo de gravedad del accidente FACTOR 3 niveles: HERIDO, MUERTO, SOLO DAÑO
DISENO Tipo de diseño en donde ocurrió el accidente Factor 13 niveles: Ciclor Ruta, Glorieta, Interseccion, Lote o Predio, Paso a nivel, Paso Elevado, Paso Inferior, Ponton, Puente, Tramo de via, Tunel, Via peatonal, Otro
DIA NOMBRE Nombre del dia de la semana del accidente Factor 7 niveles: LUNES, MARTES, MIERCOLES, JUEVES, VIERNES, SABADO, DOMINGO
MES Número de mes dentro del año del accidente integer -
CB Identificación del lote más cercano a la ubicación del accidente string Acrónimo de Comuna y barrio
COMUNA Comuna donde ocurrió el accidente Factor 21 niveles
BARRIO Barrio donde ocurrió el accidente Factor 307 niveles
ESPECIAL Fecha especial Factor niveles: SI, NO
HORA_PICO Hora pico en la ciudad (5am - 8am, 5pm-8pm) Factor niveles: SI, NO
HUMEDAD_REL Porcentaje Humedad relativa Numeric -
PLUVO Cantidad de precipitaciones [mm] Numeric -

3.3 Validación datos (Análisis descriptivo)

En los últimos años ha habido un acelerado incremento del parque automotor circulante en las ciudades del país, lo que genera problemas de movilidad, tales como la congestión y aumento de la accidentalidad. De acuerdo con Censo de Población y Vivienda 2018 realizado por el DANE en el país hay 48.258.494 habitantes de los cuales se estima el 30.3% tenía hasta Junio 20 de 2019 algún tipo de vehiculo registrado en el parque automotor donde el 17.2% de esa población tenía principalmente motocicletas, el 7.1% automoviles y el 3.3% camionetas.

En medellín, según cálculos del Registro Único Nacional de Tránsito (Runt) el parque automotor del Valle de Aburrá viene en crecimiento continuo durante los últimos años. Solo entre 2016 y 2019 el número de vehiculos se incrementó en más de 210.000 unidades, arrojando un total de 1.674.554 vehiculos en la ciudad. Desde 2014 a 2018 en la ciudad se presentaron 228.693 accidentes, en los cuales ocasionaron la muerte a 1.374 personas, 126.222 resultaron heridas y además en 101.097 accidentes dejaron solo daños materiales. El siguiente análisis descriptivo se realiza teniendo en cuenta el objetivo de este trabajo y cabe resaltar que en la preparación de los datos se eliminaron 19.948 registros repartidos en los 5 años debido a duplicados o ausencia del dato.

3.3.1 Análisis temporal

<p align = "justify> Los siguientes grárficos de barras muestran el total de accidentes para 4 resoluciones de tiempo diferentes. En general, en Medellín diariamente en promedio ocurren 114 accidentes, un promedio semanal de 787 accidentes, mensual de 3479 y al año en promedio se presentan 41.749 accidentes. Los Domingos son cuando menos ocurren accidentes y a comienzo de año también se nota una reducción de accidentes. La mayoría de accidentes diarios oscilan entre 50 y 150, sin embargo hubo días en que el numero de accidentes sobrepasó mucho o hubo pocos; el día que presentó más accidentes fue el Viernes, 4 de Agosto de 2017 con 217 accidentes y el día con menor registro fue el Viernes, 25 de Marzo de 2016 con 28 accidentes.

require(ggplot2)
require(ggpubr)

#============= Análisis Temporal

#===============
#     Dias
#===============
# Registros por día
medellin_dia <- medellin %>% group_by(DIA_NOMBRE) %>% summarize(total_accidentes = n())

# Ordenar de Lunes a Domingo
medellin_dia$DIA_NOMBRE <- ordered(medellin_dia$DIA_NOMBRE, 
                                   levels = c("LUNES","MARTES","MIERCOLES","JUEVES","VIERNES",
                                              "SABADO","DOMINGO"))


grafico_1 <- ggplot(data = medellin_dia, aes(x = DIA_NOMBRE, y = total_accidentes)) +
  geom_bar(stat = "identity", alpha = 0.65, fill = "grey20", color = "black") + 
  geom_text(aes(y = total_accidentes, label = total_accidentes),
            position = position_dodge(width = 0.8), size = 2.8, vjust = 0.5, hjust = -0.1, col = "gray10") +
  xlab("Días") + 
  ylab("Total Accidentes") + 
  ggtitle("Total Accidentes por día de la semana") +
  theme_minimal() +
  ylim(c(0,40000)) +
  coord_flip()

#=================
#     Mes
#=================
nombre_mes = c("Ene","Feb","Mar","Abr","May","Jun","Jul","Ago","Sept","Oct","Nov","Dic")
med_mes <- table(medellin$MES)
names(med_mes) <- nombre_mes
med_mes <- as.data.frame(med_mes)

grafico_2 <- ggplot(data = med_mes, aes(x = Var1, y = Freq)) + 
  geom_bar(stat = "identity", alpha = 0.7, fill = "grey20", col = "black") +
  geom_text(aes(y = Freq, label = Freq),
            position = position_dodge(width = 0.7), size = 2, vjust = -0.5, hjust = 1.5, col = "gray10") +
  xlab("Mes") +
  ylab("Total Accidentes") +
  ggtitle("Total Accidentes por mes") +
  ylim(c(0,20000)) +
  theme_minimal()


#=================
#     Anio
#=================
med_Anio <- table(medellin$PERIODO) %>% as.data.frame()

grafico_3 <- ggplot(data = med_Anio, aes(x = Var1, y = Freq)) +
  geom_bar(stat = "identity", alpha = 0.65, fill = "grey20", col = "black") + 
  geom_text(aes(y = Freq, label = Freq),
            position = position_dodge(width = 0.7), size = 3.5, vjust = 1.2, hjust = 0.5, col = "gray10") +
  xlab("Año") +
  ylab("Total Accidentes") +
  ggtitle("Total Accidentes por año") + 
  theme_minimal()

#===================
#     Semana
#===================
med_sem <- medellin %>% group_by(semana=week(medellin$FECHA)) %>% 
  summarise(Total_Accidentes = n())

grafico_4 <- ggplot(data = med_sem, aes(x = semana, y = Total_Accidentes)) +
  geom_bar(stat = "identity", alpha = 0.6, fill = "grey20", col = "black") + 
  xlab("Año") +
  ylab("Total Accidentes") +
  ggtitle("Total Accidentes por Semana") + 
  theme_minimal()

# Imprimir los 4 graficos
ggarrange(grafico_3, grafico_2, grafico_4,grafico_1, ncol = 2, nrow = 2)

#======================
# Serie temporal diaria
#=======================
med_diaria <- medellin %>% group_by(FECHA) %>% summarize(total_accidentes = n())
med_diaria$Anio <- year(med_diaria$FECHA)
med_diaria$dia <- day(med_diaria$FECHA)
med_diaria$mes <- month(med_diaria$FECHA)

plot_ly(data = med_diaria, x = ~FECHA, y = ~total_accidentes,
        type = "scatter", mode = "lines", split = ~Anio,
        line = list(width = 1)) %>%
  layout(title = 'Registros de accidentalidad diarios 2014-2018',
        xaxis = list(title = "Día"),
        yaxis = list(title = "Número de accidentes"))

El mes con menos promedio de accidentes diarios es Enero, con un promedio de 97 accidentes, mientras que Agosto presenta el mayor promedio de accidentes diarios con 121, esto teniendo en cuenta los 5 años de registros. El domingo es el día de la semana que menos accidentes se presentan con un promedio de 75 mientras que los demás días están alrededor de 100 y 140 accidentes en promedio. El día que presenta mayor variabilidad en cantidad de accidentes es el viernes, presentando varios datos atipicos por encima y por debajo.

# Gráfico de lineas Promedio diario mensual por anio
aggregate(total_accidentes~Anio*mes, data = med_diaria, FUN = mean) %>%
  plot_ly(x = ~mes,
          y = ~total_accidentes, type = "scatter", mode = "lines",
          split = ~Anio, line = list(width = 1.5)) %>%
  layout(title = 'Promedio diario mensual por año',
         xaxis = list(title = "Mes"),
         yaxis = list(title = "Número de accidentes"))
# Boxplot por dia de la semana
med_diaria$semana <- weekdays(med_diaria$FECHA)
med_diaria$semana <- factor(med_diaria$semana,
                             levels = c("lunes","martes","miércoles","jueves","viernes","sábado","domingo"))

plot_ly(data = med_diaria, x = ~semana, y = ~total_accidentes,
        type = "box", marker = list(color = 'grey24')) %>%
  layout(title = 'Accidentes por día de la semana',
         xaxis = list(title = "Día de la semana"),
         yaxis = list(title = "Accidentes"))

3.3.2 Análisis Geoespacial Barrios y Comunas

Se analiza la variable barrio y comuna, La Candelaria es la Comuna que mas presenta accidentes entre 2014 y 2018 seguida por Poblado. El Barrio donde más se presenta accidentes es justamente es el barrio La Candelaria, también se puede notar del mapa que la zona centro es la más afectada por los accidentes y las vias principales.

# Accidentes por comuna
medellin_comuna <- medellin %>% group_by(COMUNA) %>%summarize(accidentes = n())

ggplot(data = medellin_comuna, aes(x = reorder(COMUNA,+accidentes), y = accidentes)) +
  geom_bar(stat = "identity", position = "dodge", fill = "grey24", color = "black", alpha = 0.6) +
  geom_text(aes(y = accidentes, label = accidentes),
                position = position_dodge(width = 0.7), size = 3.5, vjust = 0.5, hjust = -0.1, col = "black") +
  xlab("Comuna") + 
  ylab("Total accidentes") +
  ggtitle("Total accidentes por comuna 2014-2018") +
  ylim(c(0,50000)) +
  theme_minimal() +
  coord_flip()

#================================================
#                     MAPA
#================================================
require(leaflet)
# Mapa para todos los barrios
nueva_base <- medellin %>% filter(PERIODO >= 2014 & PERIODO <= 2018) %>% 
  group_by(CB) %>%
  summarise(accidentes = n()) %>%
  ungroup()

nuevo_mapa <- inner_join(catastro, nueva_base, by = c("CODIGO" = "CB"))

# Paleta de color
mypal <- colorNumeric(palette = c("#000000","#280100","#3D0201","#630201","#890100","#B00100","#DD0100","#F50201",
                                   "#FF5F5E","#FF7A79","#FF9796","#FEB1B0","#FDC9C8", "#FFE5E4"), domain = nuevo_mapa$accidentes, reverse = T)

# Crear mapa
leaflet() %>% addPolygons(data = nuevo_mapa, color = "#0A0A0A", opacity = 0.6, weight = 1, fillColor = ~mypal(nuevo_mapa$accidentes),
                          fillOpacity = 0.6, label = ~NOMBRE_BAR,
                          highlightOptions = highlightOptions(color = "black", weight = 3, bringToFront = T, opacity = 1),
                          popup = paste("Barrio: ", nuevo_mapa$NOMBRE_BAR, "<br>", "Accidentes: ", nuevo_mapa$accidentes, "<br>")) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addLegend(position = "bottomright", pal = mypal, values = nuevo_mapa$accidentes, title = "Accidentes", opacity = 0.6)

3.3.3 Análisis de factores de riesgo

Solo el 0.6% de los accidentes ha dejado muertos el otro 99.4% ocasionan Heridos o solo daños. El 88% de los accidentes ocurren en un Trammo de vía o intersección y la clase de accidente más frecuente es el Choque, además los choques son los que generan más daños y el atropello es el que genera más muertos junto a los choques, la frecuencia de heridos se nota en todos los tipos de accidentes.

# ------------ Gravedad del accidente
med_gravedad <- table(medellin$GRAVEDAD) %>% prop.table() %>% as.data.frame
med_gravedad$Var1 <- ordered(med_gravedad$Var1, levels = c("HERIDO","SOLO DAÑOS","MUERTO"))

grafico_5 <- ggplot(data = med_gravedad, aes(x = Var1, y = Freq)) +
  geom_bar(stat = "identity", alpha = 0.7, fill = "firebrick", color = "black") + 
  geom_text(aes(y = Freq, label = paste(round(Freq*100,1),"%")),
            position = position_dodge(width = 0.8), size = 4, vjust = -0.5, hjust = 0.5, col = "gray10") +
  xlab("Gravedad") + 
  ylab("Frecuencia (%)") + 
  ggtitle("Total Accidentes por gravedad") +
  theme_minimal()

#-------- Diseño de la via en el accidente
med_diseno <- table(medellin$DISENO) %>% prop.table() %>% as.data.frame
med_diseno$Var1 <- ordered(med_diseno$Var1, levels = c("Via peatonal","Tunel","Ponton","Paso a Nivel","Paso Inferior","Puente","Ciclo Ruta",
                                                       "Paso Elevado","Otro","Glorieta","Lote o Predio","Interseccion","Tramo de via"))

grafico_6 <- ggplot(data = med_diseno, aes(x = Var1, y = Freq)) +
  geom_bar(stat = "identity", alpha = 0.7, fill = "firebrick", color = "black") + 
  geom_text(aes(y = Freq, label = paste(round(Freq*100,2),"%")),
            position = position_dodge(width = 0.8), size = 3.2, vjust = 0.5, hjust = -0.2, col = "gray10") +
  xlab("Diseño") + 
  ylab("Frecuencia (%)") + 
  ggtitle("Total Accidentes por diseño de la via") +
  theme_minimal() + 
  coord_flip() +
  ylim(0,1)

#------------ Clase Accidente
med_clase <- table(medellin$CLASE) %>% prop.table() %>% as.data.frame
med_clase$Var1 <- ordered(med_clase$Var1, levels = c("Choque","Otro","Atropello","Caida Ocupante","Volcamiento","Incendio","Choque y Atropello"))

grafico_7 <- ggplot(data = med_clase, aes(x = Var1, y = Freq)) +
  geom_bar(stat = "identity", alpha = 0.7, fill = "firebrick", color = "black") + 
  geom_text(aes(y = Freq, label = paste(round(Freq*100,3),"%")),
            position = position_dodge(width = 0.8), size = 4, vjust = -0.5, hjust = 0.5, col = "gray10") +
  xlab("Clase") + 
  ylab("Frecuencia relativa (Accidentes)") + 
  ggtitle("Total Accidentes por Clase") +
  theme_minimal() + 
  ylim(0,0.8)

# Imprimir los 3 graficos
ggarrange(grafico_5, grafico_6, ncol = 2, nrow = 1)

grafico_7

# Clase vs gravedad
med_clase_gravedad <- table(medellin$CLASE, medellin$GRAVEDAD) %>% prop.table() %>% as.data.frame()
med_clase_gravedad$Var1 <- ordered(med_clase_gravedad$Var1, levels = c("Choque y Atropello","Incendio","Volcamiento","Caida Ocupante","Atropello","Otro","Choque"))
names(med_clase_gravedad) <- c("Var1","Gravedad","Freq")

ggplot(data = med_clase_gravedad, aes(x = Var1, y = Freq, fill = Gravedad)) + 
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(y = Freq, label = paste(round(Freq*100,1),"%")),
            position = position_dodge(width = 0.8), size = 3, vjust = 0.5, hjust = -0.2, col = "gray10") +
  theme_minimal() +
  ylim(0,0.5) +
  xlab("Frecuencia (%)") +
  ylab("Clase de Accidente") +
  coord_flip()

3.4 Agrupamiento de barrios

Con el fin de agrupar barrios según sus caracteristicas de accidentalidad se usa k-means para buscar los barrios mas cercanos a dichas caracteristicas. Para este trabajo las variables que caracterizan a los barrios son las siguientes:

  • Total accidentes por barrio
  • Total de muertes en accidentes por barrio
  • Total de accidentes en fechas especiales
  • Total de accidentes en hora pico

Para obtener el número k optimo que reduce la distancia entre cluster, se realiza un gráfico de codos en el que se nota que el número k optimo es 4 u 5, asi que elegimos un \(k = 5\) dado la cantidad de barrios de la ciudad.

# Seleccionar las variables que se usan
agrupamiento <- medellin[,c("FECHA","BARRIO","COMUNA","PERIODO","CB","ESPECIAL","HORA_PICO",
                            "GRAVEDAD")]

# Agrupar las variables por barrijo
data_agrupamiento <- agrupamiento %>% filter(PERIODO >= 2014 & PERIODO <= 2018) %>% 
  group_by(CB) %>%
  summarise(total_accidentes = n(), 
            Acc_Especiales = sum(ESPECIAL == "SI"),
            Acc_Hpico = sum(HORA_PICO == "SI"),
            Acc_Muerto = sum(GRAVEDAD == "MUERTO")) %>%
  ungroup()

# Grafico de codos para k optimo
wssplot <- function(data, nc = 15, seed = 1234){
  wss <- (nrow(data)-1)*sum(apply(data,2,var))
  for (i in 2:nc){
    set.seed(seed)
    wss[i]<- sum(kmeans(data, centers = i)$withinss)}
  plot(1:nc, wss, type = "b", xlab = "Number of Clusters",
       ylab = "Within groups sum of squares", main = "Busqueda de numero de cluster optimo")}
wssplot(data_agrupamiento[,-1], nc = 10)

Ahora, teniendo en cuenta que el numero de grupos o cluster optimo es 5, entonces se procede a usar kmeans de R para agrupar los barrios en 5 grupos y se realiza un grafico con sus dos primeras componentes principales para visualizar el agrupamiento. En la gráfica se puede notar que la dimension 1 que corresponde a el total de accidentes por barrio explica la variabilidad en un 79.5%.. Los datos de agrupamiento se guardan en un archivo csv llamado Agrupamiento_barrios.csv

require(factoextra)
require(tidyverse)
require(ggplot2)

# Ejecucion modelo para agrupamiento con kmenas
set.seed(1017239703)
barrios <- data_agrupamiento$CB
k.means.fit <- kmeans(data_agrupamiento[,-1], 5)
d <- k.means.fit$cluster
data_agrupamiento$Grupos <- d

# Visualizacion de grupos con sus dos primeras componentes principales
agrup_grafica <- data_agrupamiento[,-1]
row.names(agrup_grafica) <- barrios
set.seed(123456789)
fviz_cluster(object = k.means.fit, data = agrup_grafica,
             ellipse.type = "convex", repel = F, show.clust.cent = FALSE,
             labelsize = 8) +
  labs(title = "Clasificación Jerárquica Barrios Medellín",
       subtitle = "K Nearest Neighbor") +
  theme_bw() +
  theme(legend.position = "bottom")

Ahora, para visualizar la agrupacion en el mapa haremos uso del shape file de catastro y de leaflet para dibujar el mapa y con colores representar los grupos seleccionados. En el mapa se puede notar claramente como los barrios del centro son los mas afectados por la accidentalidad y los barrios cercanos a las vias rapidas de la ciudad también son muy afectados. La idea del mapa es poder servir como gestor para implementación de estrategias que puedan reducir la accidentalidad de una forma más precisa, como ejemplo en la comuna La Candelaria que es la mas afectada por los accidentes y la mortalidad pero que cuenta con 20 barrios en los que se encuentra el Centro administrativo o el barrio Las Palmes que hacen parte del grupo de los mas bajos accidentes, pero se encuentra, Guayaquil, La Candelaria, Barrio Colón, Perpetuo Socorro y San Benito que hacen parte del selecto grupo de los de mayor accidentalidad en Medellín.

setwd("Datos")
data_agrupamiento <- read.csv(file = "Otros/Agrupamiento_barrios.csv", sep = ",", encoding = "UTF-8", colClasses = "character")
# Mapa para todos los barrios
dato_mapa <- inner_join(catastro, data_agrupamiento, by = c("CODIGO" = "CB"))


# Definicion de colores
thalia <- c("#FF0000", "#CCFF00", "#00FF66", "#0066FF", "#CC00FF")
dato_mapa$colores <- ifelse(dato_mapa$Grupos == "Grupo 1", "#FF0000",
                            ifelse(dato_mapa$Grupos == "Grupo 2", "#CCFF00",
                                   ifelse(dato_mapa$Grupos == "Grupo 3", "#00FF66",
                                          ifelse(dato_mapa$Grupos == "Grupo 4", "#0066FF",
                                                 ifelse(dato_mapa$Grupos == "Grupo 5","#CC00FF",0)))))



# Crear mapa
leaflet() %>% addPolygons(data = dato_mapa, opacity = 0.4, color = "#545454",weight = 1, fillColor = dato_mapa$colores,
                          fillOpacity = 0.4, label = ~NOMBRE_BAR,
                          highlightOptions = highlightOptions(color = "#262626", weight = 3, bringToFront = T, opacity = 1),
                          popup = paste("Barrio: ", dato_mapa$NOMBRE_BAR, "<br>", "Grupo: ", dato_mapa$Grupos)) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addLegend(position = "bottomright", colors = thalia, labels = c("Grupo 1: Muy alta accidentalidad y mortalidad", "Grupo 2: Poca accidentalidad", "Grupo 3: Moderada accidentalidad", "Grupo 4: Alta accidentalidad", "Grupo 5: Muy poca accidentalidad"))

3.5 Prototipo, Validación y evaluación

Para Predecir la accidentalidad de la ciudad de medellín se ha dispuesto de la informacion suministrada por el IDEAM (Instituto de Hidrología, Meteorología y Estudios Ambientales), con esta se recaudo la medida de la pluviometría diaria y la humedad relativa a nivel del suelo en la ciudad de medellin.

A partir de esta información y la obtenida por el canal de datos abiertos de movilidad que publica la Alcaldía de Medellín en el portal GeoMedellín donde se exhiben los diferentes accidentes diarios y medidas del mismo,junto con la compilacion de las fechas festivas colombianas y generando un atributo para determinar si el accidente ocurrio en medio de horas pico, teniendo en cuenta que se desea predecir la accidentalidad por tipo de accidente a nivel semanal, mensual y diario, se procedio a separar esta base de datos en dos, de la historia de accidentalidad que comprende desde el 2014 hasta el 2017, y la del 2018, esta ultima para realizar validaciones de los modelos planteados, esta base de datos la llamamos accidentalidad_medellin_final.csv y puede descargarla de github.

Asi, se procedio a agrupar todos los accidentes sucedidos en un mismo Dia, o en una misma semana, o en un mismo mes,junto con el tipo de accidente (Choque, Atropello, Caida Ocupante, Volcamiento, Incendio, Otro) usando funciones de agregacion para los diferentes atributos expuestos anteriormente, puede ver el script y las bases de datos en aqui.

Ya agrupadas y listas las bases de datos por mes,semana y dia al igual que por tipo de accidentes, procederemos a explicar brevemente los modelos que se utilizaron para generar las predicciones del año 2018 y sus desempeños.

Modelación

A continuación se dará una breve descripción de la manera en que se trataron y usaron las bases de datos para entrenar modelos predictivo sobre la cantidad de accidentes en un periodo de tiempo, para esto se evaluaron 5 tipos de modelos:

  • XGBoost
  • Regresión lineal generalizada con distribución exponencial usando GAMLSS
  • Random Forest
  • K-Nearest Neighbor Regression
  • Redes Neuronales

También es importante tener en cuenta las resoluciones temporales que se trabajarán los modelos y son las siguientes:

  • Diario
  • Semanal
  • Mensual

Los registros faltantes en las variables del clima se les realizó imputación basado en su media para la Humedad relativa y valor de cero para la Pluviometria. A continuación se presenta como queda conformado el conjunto de entrenamiento para la resolución semanal, de esta misma forma queda construido el nivel diario y mensual para los conjuntos de entrenamiento y validación.

require(lubridate)
require(tidyverse)
require(fastDummies)
require(xgboost)
require(Matrix)
require(gamlss)
require(skimr)
require(ranger)
library(FNN)
require(neuralnet)
require(gotop)
use_gotop()
#==============================
#  Imputación datos faltantes
#===============================
medellin$HUMEDAD_REL[is.na(medellin$HUMEDAD_REL)] <- mean(medellin$HUMEDAD_REL, na.rm = T)
medellin$PLUVIO[is.na(medellin$PLUVIO)] <- 0

#=======================================================
#   Elección conjunto de entrenamiento y validación
#=======================================================
Train <- medellin[medellin$PERIODO %in% c(2014,2015,2016,2017),]
Test <- medellin[medellin$PERIODO == 2018,]

#=========================
#    Puesta a Punto
#==========================

# Convertir ESPECIAL y HORA PICO en binarias
Train$ESPECIAL <- ifelse(Train$ESPECIAL=='SI',1,0)
Train$HORA_PICO <- ifelse(Train$HORA_PICO=='SI',1,0)

Test$ESPECIAL <- ifelse(Test$ESPECIAL=='SI',1,0)
Test$HORA_PICO <- ifelse(Test$HORA_PICO=='SI',1,0)

Train$CLASE <- as.factor(Train$CLASE)
levels(Train$CLASE) <- c('Atropello','Caida_ocupante', 'Choque','Choque',
                         'Incendio','Otro','Volcamiento')
Test$CLASE <- as.factor(Test$CLASE)
levels(Test$CLASE) <- c('Atropello','Caida_ocupante', 'Choque','Choque',
                        'Incendio','Otro','Volcamiento')

#====================================================
#  Agrupamiento variables temporales para los modelos
#====================================================
# Grupo1: Mensual   
Mensual <- Train %>% group_by(Mes = month(FECHA),CLASE) %>% 
  summarise(cantidad = length(FECHA),
            festivos = sum(ESPECIAL),
            Hora_Pico = sum(HORA_PICO),
            ValorHumR = mean(HUMEDAD_REL, na.rm = T),
            Pluvial  = max(PLUVIO)) %>% 
  dummy_cols(select_columns = c('CLASE'),
             remove_selected_columns = T)

MensualTs <- Test %>% group_by(Mes = month(FECHA),CLASE) %>% 
  summarise(cantidad = length(FECHA),
            festivos = sum(ESPECIAL),
            Hora_Pico = sum(HORA_PICO),
            ValorHumR = mean(HUMEDAD_REL, na.rm = T),
            Pluvial  = max(PLUVIO)) %>% 
  dummy_cols(select_columns = c('CLASE'),
             remove_selected_columns = T)   

# Grupo 2: Diario
Diario <- Train %>% group_by(Dia = day(FECHA),CLASE) %>% 
  summarise(cantidad = length(FECHA),
            festivos = sum(ESPECIAL),
            Hora_Pico = sum(HORA_PICO),
            ValorHumR = mean(HUMEDAD_REL, na.rm = T),
            Pluvial  = max(PLUVIO)) %>% 
  dummy_cols(select_columns = c('CLASE'),
             remove_selected_columns = T)

DiarioTs <- Test %>% group_by(Dia = day(FECHA),CLASE) %>% 
  summarise(cantidad = length(FECHA),
            festivos = sum(ESPECIAL),
            Hora_Pico = sum(HORA_PICO),
            ValorHumR = mean(HUMEDAD_REL, na.rm = T),
            Pluvial  = max(PLUVIO)) %>% 
  dummy_cols(select_columns = c('CLASE'),
             remove_selected_columns = T)

# Grupo 3: Semanal
Semanal <- Train %>% group_by(Semana = week(FECHA),CLASE) %>% 
  summarise(cantidad = length(FECHA),
            festivos = sum(ESPECIAL),
            Hora_Pico = sum(HORA_PICO),
            ValorHumR = mean(HUMEDAD_REL, na.rm = T),
            Pluvial  = max(PLUVIO)) %>% 
  dummy_cols(select_columns = c('CLASE'),
             remove_selected_columns = T)

SemanalTs <- Test %>% group_by(Semana = week(FECHA),CLASE) %>% 
  summarise(cantidad = length(FECHA),
            festivos = sum(ESPECIAL),
            Hora_Pico = sum(HORA_PICO),
            ValorHumR = mean(HUMEDAD_REL, na.rm = T),
            Pluvial  = max(PLUVIO)) %>% 
  dummy_cols(select_columns = c('CLASE'),
             remove_selected_columns = T)

head(Semanal)

Ahora, se describe el proceso uno a uno de los modelos utilizados y su implemantación en los datos:

3.5.1 XGBoost

XGBoost es una biblioteca optimizada de aumento de gradiente distribuida diseñada para ser altamente eficiente, flexible y portátil. Implementa algoritmos de aprendizaje automático en el marco de Gradient Boosting.

Boosting es una metodologia para mejorar la prediccion proporcionada por un arbol de decision. Como Bagging, Boosting es un procedimiento general que puede ser aplicado a muchos metodos de aprendizaje estadistico para regresion y clasicacion. En Boosting en vez de ajustar un arbol de decision grande, lo cual tiene el potencial de de sobre ajustar los datos, en forma secuencial se adicionan arboles donde cada uno es ajustado a los residuos del anterior.

Tome \(\hat{f}(x)=0\) y \(r_i=y_i\) para todo i en el conjunto de entrenamiento.

ahora, para un b=1,2,3,…,B se repiten los siguientes pasos:

  • Ajuste un arbol con d bifurcaciones (d+1 nodos terminales) a los datos (\(X,r\)).

  • Actualice \(\hat{f}\) adicionando una version ‘contraida’ del arbol nuevo, es decir:

\[\hat{f} \rightarrow \hat{f}+\lambda\hat{f^b}(x)\] - Actualice los residuos:

\[r_i=r_i-\lambda\hat{f^b}(x_i)\] Finalmente se obtiene un modelo Boosted

\[\hat{f}=\sum_{b=1}^B\lambda\hat{f^b}(x)\]

Ahora, se genero un modelo XGboost que reducia la funcion de costo de la raiz del MSE, con un parametro de aprendizaje de \(\lambda=1\) y 100 arboles.

La forma en la que se entreno el modelo en R es como sigue:

#============================
#          XGBoost
#============================
param1 <- list(booster = "gblinear", nthread = 100, eta = 0.1,
               gamma = 0, max = 10, min_child_weight = 1,
               max_delta_step = 2, lambda = 1,
               alpha = 0, three_method = "exact",
               objective = "reg:squarederror",
               eval_metric = "rmse")  # Parametros del Boosting

#-----------------
# XGboost: Mensual
#-----------------
dtrain <- xgb.DMatrix(as.matrix(Mensual[,-c(2)]),label = Mensual$cantidad)
dtest <- xgb.DMatrix(as.matrix(MensualTs[,-c(2)]),label = MensualTs$cantidad)

watchlist <- list(train=dtrain,test=dtest)

ModMesXGboost <- xgb.train(data = dtrain,nrounds = 200,
                           watchlist = watchlist,params = param1,
                           max.depth=2)
#------------------
# XGboost: Semanal
#------------------
dtrain <- xgb.DMatrix(as.matrix(Semanal[,-c(2)]),label = Semanal$cantidad)
dtest <- xgb.DMatrix(as.matrix(SemanalTs[,-c(2)]),label = SemanalTs$cantidad)

watchlist <- list(train=dtrain,test=dtest)

ModSemXGboost <- xgb.train(data = dtrain,nrounds = 100,
                           watchlist = watchlist,params = param1,
                           max.depth=2)
#-----------------
# XGboost: Diario
#-----------------
dtrain <- xgb.DMatrix(as.matrix(Diario[,-c(2)]),label = Diario$cantidad)
dtest <- xgb.DMatrix(as.matrix(DiarioTs[,-c(2)]),label = DiarioTs$cantidad)

watchlist <- list(train=dtrain,test=dtest)

ModDiaXGboost <- xgb.train(data = dtrain,nrounds = 100,
                           watchlist = watchlist,params = param1,
                           max.depth=2)

#------------------------
# Validacion MSE XGboost
#------------------------
require(Matrix)

MSEXGboostMes <- mean((predict(ModMesXGboost,
                               as(as.matrix(MensualTs[,-2]),"dgCMatrix"))-MensualTs$cantidad)^2)

MSEXGboostSem <- mean((predict(ModSemXGboost,
                               as(as.matrix(SemanalTs[,-2]),"dgCMatrix"))-SemanalTs$cantidad)^2)

MSEXGboostDia <- mean((predict(ModDiaXGboost,
                               as(as.matrix(DiarioTs[,-2]),"dgCMatrix"))-DiarioTs$cantidad)^2)

finalmente se entreno el modelo para las bases de agrupadas por Mes, Semana y Dia asi como por tipo de accidente, se obtuvieron los siguientes resultados en la validacion cruzada por el metodo de conjunto de validacion sobre el año 2018.

  • MSE Periodo Mensual: 7882.709
  • MSE Periodo Semanal: 507.5032
  • MSE Periodo Diario: 1522.026

3.5.2 Regresión lineal generalizada usando GAMLSS

Los GAMLSS (Generalized Additive Models for Location,Scale and Shape) son modelos de regresión distributiva univariante, donde todos los parámetros de la distribución asumida para la respuesta se pueden modelar como funciones aditivas de las variables explicativas.

Es una mejora del modelo GAM (Generalized Additive Model), definido asi:

\(Y_i = \beta_0+f_1(X_{i1})+f_2(X_{i2})+...+f_p(X_{ip})+\epsilon\)

Para el modelo GAMLSS consideramos los siguientes vectores de parametros:

  • \(\mu :\) Parámetro relacionados con la ubicación, por ejemplo, media
  • \(\sigma :\) Parámetro relacionados con la escala, por ejemplo, varianza
  • \(\nu :\) Parámetro relacionados con la forma, por ejemplo, asimetría
  • \(\tau :\) Parámetro relacionados con la forma, por ejemplo, curtosis

Asumiendo que Y es la variable respuesta y \(D\) es alguna distribucion, el modelo GAMlSS se plantea como:

\[Y \sim D(\mu,\sigma,\nu,\tau)\] \[g_1(\mu) = \beta_{01}+f_{11}(X_1)+f_{12}(X_2)+ ... + f_{1p}(X_p)\] \[g_2(\mu) = \beta_{02}+f_{21}(X_1)+f_{22}(X_2)+ ... + f_{2p}(X_p)\] \[g_3(\mu) = \beta_{03}+f_{31}(X_1)+f_{32}(X_2)+ ... + f_{3p}(X_p)\] \[g_4(\mu) = \beta_{04}+f_{41}(X_1)+f_{42}(X_2)+ ... + f_{4p}(X_p)\]

Con lo que a diferentes funciones \(f\) o tambien llamadas links y una distribucion \(D\) podremos generar un modelo que prediga la cantidad de accidentes en un tiempo determinado (experimento \(Poisson\))

Para nuestro caso particular utilizamos como funcion \(D\) una \(Poisson\) y el paquete gamlss posee una funcion de ajuste aproximado para \(D\) fitDist() con la que tambien se modelo una \(Exponential\), siendo esta ultima la que mejor desempeño en prueba obtuvo con los siguientes valores por el metodo de validacion en el año 2018.

La forma en la que se entreno el modelo en R es como sigue:

#=========================================================================
#               Regresión lineal generalizada usando GAMLSS
#==========================================================================

#------------------
# RLGamlss: Mensual
#-------------------
ModMenGam <- gamlss(formula = cantidad~.,
                    sigma.formula = ~.,
                    data = Mensual,
                    family = EXP)
#-------------------
# RLGamlss: Semanal
#--------------------
ModSemGam <- gamlss(formula = cantidad~.,
                    sigma.formula = ~.,
                    data = na.omit(Semanal),
                    family = EXP)
#--------------------
# RLGamlss: Diario
#--------------------
ModDiaGam <- gamlss(formula = cantidad~.,
                    sigma.formula = ~.,
                    data = na.omit(Diario),
                    family = EXP)

#-------------------------------
# Validacion MSE para RLGamlss
#-------------------------------
MSEgamlssMes <- mean((predictAll(ModMenGam, MensualTs)$mu-MensualTs$cantidad)^2)

MSEgamlssSem <- mean((predictAll(ModSemGam, SemanalTs)$mu-SemanalTs$cantidad)^2)

MSEgamlssDia <- mean((predictAll(ModDiaGam, DiarioTs)$mu-DiarioTs$cantidad)^2)
  • MSE Periodo Mensual: 1151768
  • MSE Periodo Semanal: 22477.15
  • MSE Periodo Diario: 96030.66

3.5.3 Random Forest

Es un metodo de ensamble de arboles de clasificacion o regresion, donde se consideran arboles sin poda cada arbol es entrenado sobre una muestra boostrap de la base de interes, y a su vez en cada split del arbol se utiliza una muestra aleatoria de los atributos en la base de interes.

En nuestro modelamiento se utilizaron 350 arboles de regresion y muestras de 4 variables en cada split de la cantidad de accidentes en cierto periodo respecto a las demas covariables, ademas se opto por una libreria considerablemente mas optimizada ranger lo que permite realizar muchas pruebas de bosques con diferentes valores para los parametros en poco tiempo.

asi, nuestro bosque aleatorio, fue entrenado en base a datos del 2014 al 2017 y validado con datos del año 2018, obteniendo las metricas del MSE por conjunto de validacion del año 2018 como sigue.

La forma en la que se entreno el modelo en R es como sigue:

#====================================
#         Random Forest
#====================================

#------------------------
# Random Forest: Mensual
#------------------------
ModMenRanger <- ranger(cantidad~., 
                       data = na.omit(Mensual),
                       num.trees = 350, mtry = 4)

predRFMes <- predict(ModMenRanger, na.omit(MensualTs))


#-----------------------
# Random Forest: Semanal
#-----------------------
ModMenRanger <- ranger(cantidad~., 
                       data = na.omit(Semanal),
                       num.trees = 350,mtry = 4)

predRFSem <- predict(ModMenRanger, na.omit(SemanalTs))

#-----------------------
# Random Forest: Diario
#-----------------------
ModMenRanger <- ranger(cantidad~., 
                       data = na.omit(Diario),
                       num.trees = 350,mtry = 4)

predRFDia <- predict(ModMenRanger, DiarioTs)


#---------------------------------
#  Validacion MSE Random Forest
#---------------------------------
MSERangerMes <- mean(((predRFMes$predictions)-MensualTs$cantidad)^2)

MSERangerSem <- mean(((predRFSem$predictions)-SemanalTs$cantidad)^2)

MSERangerDia <- mean(((predRFDia$predictions)-DiarioTs$cantidad)^2)
  • MSE Periodo Mensual : 1404420
  • MSE Periodo Semanal : 35593.04
  • MSE Periodo Diario : 301954.1

3.5.4 K-Nearest Neighbor Regression

Este metodo permite obtener el valor de una prediccion a travez de la media de sus k vecinos mas cercanos (principalmente metricas euclidianas), es decir, para hallar el valor predicho de una observacion \(X_i\) en \(y_i\) y dado un valor de k, se promedian los valores \(y\) de las k observaciones mas cercanas a la de interes, este promedio es el valor predicho para esta observacion.

En nuestro caso se probaron valores de k = 1,2,3,…,50 para cada base de datos agrupada, asi se encontro el K-optimo y con este se obtuvo los siguientes resultados de validacion por el metodo de conjunto de validacion con datos del 2018.

La forma en la que se entreno el modelo en R es como sigue:

#========================================
#       K-Nearest Neighbor Regression
#========================================

#-----------------
#   KNN: Mensual
#-----------------
MSEKnnMes <- vector()
for(i in 1:50){
  ModMesKnn <- knn.reg(train = na.omit(Mensual[,-c(2,3)]),
                       test = na.omit(MensualTs[,-c(2,3)]),
                       y = Mensual$cantidad,k=i)
  
  MSEKnnMes[i] <- mean(((ModMesKnn$pred)-na.omit(MensualTs)$cantidad)^2)
}
plot(1:50,MSEKnnMes, type = 'l')
paste("el valor de K que minimisa el MSE es:",which.min(MSEKnnMes))

#-------------
# KNN: Semanal
#-------------
MSEKnnSem <- vector()
for(i in 1:50){
  ModSemKnn <- knn.reg(train = na.omit(Semanal[,-c(2,3)]),
                       test = na.omit(SemanalTs[,-c(2,3)]),
                       y = Semanal$cantidad,k=i)
  
  MSEKnnSem[i] <- mean(((ModSemKnn$pred)-na.omit(SemanalTs)$cantidad)^2)
}
plot(1:50,MSEKnnSem, type = 'l')
paste("el valor de K que minimisa el MSE es:",which.min(MSEKnnSem))

#-------------
# KNN: Diario
#-------------
MSEKnnDia <- vector()
for(i in 1:50){
  ModDiaKnn <- knn.reg(train = na.omit(Diario[,-c(2,3)]),
                       test = na.omit(DiarioTs[,-c(2,3)]),
                       y = Diario$cantidad,k=i)
  
  MSEKnnDia[i] <- mean(((ModDiaKnn$pred)-na.omit(DiarioTs)$cantidad)^2)
}
plot(1:50,MSEKnnDia, type = 'l')
paste("el valor de K que minimisa el MSE es:",which.min(MSEKnnDia))
  • MSE Periodo Mensual: 98177.22
  • MSE Periodo Semanal: 4888.43
  • MSE Periodo Diario: 16135.33

3.5.5 Redes Neuronales

Las redes neuronales son un modelo computacional vagamente inspirado en el comportamiento observado en su homólogo biológico. Consiste en un conjunto de unidades, llamadas neuronas artificiales, conectadas entre sí para transmitirse señales. La información de entrada atraviesa la red neuronal (donde se somete a diversas operaciones) produciendo unos valores de salida.

El algoritmo basicamente genera en principio unos pesos \(W_1\) aleatorios entre las capas definidas por el usuario. De forma resumida, estos pesos se multiplican por cada registro y este resultado se le aplica una funcion no lineal a la que llamaremos activacion, y con esto obtendremos los valores de la primer capa de la red, a estos valores se le asignan unos pesos \(W_2\) aleatorios para realizar de nuevo la operacion anteriormente descrita y obtener la segunda capa de la red asi sucesivamente entre todas las capas de la estructura definida por el usuario.

Cuando se llega a la ultima capa llamada capa de salida(definida por el usuario) y despues de haber calculado todas las neuronas y sus activaciones (llamado propagacion hacia adelante), el algoritmo calcula una funcion de costo o funcion objetivo, con la intencion de reducir esta funcion (generalmente una distancia ponderada de el valor de la capa de salida a un vector \(y\) llamado Target) se recurre a recalcular los \(W\) para que se minimice la funcion de costo, a esto se le llama propagacion hacia atras, cabe aclarar que en cada activacion de las neuronas se generan vectores \(b\) que se suman a cada neurona.

En nuestro caso particular y notando que el metodo de redes es computacionalmente demandante, generamos las redes neuronales que se podian correr en nuestra maquina, es decir utilizamos una arquitectura de red tal que pudiese correrse en nuestra maquina y fueron las siguientes:

Mensual
Red neuronal Mensual
ModMenNN <- neuralnet(cantidad~.,data = na.omit(Mensual),
              hidden = c(3,2))
Semanal
Red neuronal Semanal
ModSemNN <- neuralnet(cantidad~.,data = na.omit(Semanal),
              hidden = c(3,2))
Diario
Red neuronal Diaria
ModDiaNN <- neuralnet(cantidad~.,data = na.omit(Diario),
              hidden = c(2),stepmax = 5000000)

Con estas calculamos los MSE de validacion para el periodo 2018, obteniendo:

  • MSE Periodo Mensual: 11520464
  • MSE Periodo Semanal: No pudo ser entrenado
  • MSE Periodo Diario: 1356718

3.5.5 Resumen desempeño modelos

A continuación verá un resumen gráfico de las validaciones de los modelos teniendo como medida la Raiz del MSE. Se resalta la notoriedad de los modelos KNN y XGboost en la reducción del MSE sobre todos los demás usando como validación el 2018, siendo el XGboost el que logra mejor desempeño para los 3 periodos, por lo tanto este será el modelo para realizar las predicciones.

Modelo <- c("XGBoost","XGBoost","XGBoost",
            "GAMLSS","GAMLSS","GAMLSS",
            "RandomForest","RandomForest","RandomForest",
            "KNN","KNN","KNN",
            "RedNeuronal","RedNeuronal","RedNeuronal")
Periodo <- c("Mensual","Semanal","Diario",
             "Mensual","Semanal","Diario",
             "Mensual","Semanal","Diario",
             "Mensual","Semanal","Diario",
             "Mensual","Semanal","Diario")
MSE <- c(7882.709,507.5032,1522.026,
         1151768,22477.15,96030.66,
         1404420,35593.04,301954.1,
         98177.22,4888.43,16135.33,
         11520464,0,1356718)

Evaluacion <- data.frame(Modelo = Modelo, Periodo = Periodo, MSE = MSE)

ggplot(data = Evaluacion, aes(x = Modelo, y = sqrt(MSE), fill = Periodo)) +
  geom_bar(stat = "identity", position = "dodge") +
  theme_minimal() +
  xlab("Modelo") +
  ylab("RMSE")

3.5.6 Modelo generalizado para predicciones con fechas dia a dia

Dado que en los anteriores modelos se obtiene con los accidentes agrupados por dias, semana y mes, se necesita un modelo que relacione una a uno los dias teniendo en cuenta las fechas. Dado que los accidentes representa un conteo, se entrenará un modelo de regresión generalizada usando la distribución Poisson, para las variables FECHA, CLASE, DIA_NOMBRE, PERIODO, ESPECIAL.

#=============================
# MODELO GENERALIZADO POISSON
#=============================
Train$CLASE <- toupper(Train$CLASE)
Test$CLASE <- toupper(Test$CLASE)

grupo_1 <- Train %>% group_by(FECHA, CLASE, DIA_NOMBRE, PERIODO, ESPECIAL) %>%
  summarise(Accidentes = n()) %>%
  arrange(FECHA)


# Modelo Poisson
modelo_poisson <- glm(Accidentes~FECHA+CLASE+DIA_NOMBRE+PERIODO+ESPECIAL,
                      data = grupo_1, family = "poisson")
summary(modelo_poisson)
## 
## Call:
## glm(formula = Accidentes ~ FECHA + CLASE + DIA_NOMBRE + PERIODO + 
##     ESPECIAL, family = "poisson", data = grupo_1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -8.0990  -1.0024  -0.0941   0.7931   6.0641  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          1.199e+02  1.722e+01   6.960 3.39e-12 ***
## FECHA                1.877e-04  2.315e-05   8.110 5.06e-16 ***
## CLASECAIDA_OCUPANTE -1.345e-01  1.122e-02 -11.987  < 2e-16 ***
## CLASECHOQUE          1.890e+00  8.221e-03 229.948  < 2e-16 ***
## CLASEINCENDIO       -2.465e+00  2.427e-01 -10.159  < 2e-16 ***
## CLASEOTRO            6.831e-02  1.066e-02   6.411 1.45e-10 ***
## CLASEVOLCAMIENTO    -1.093e+00  1.558e-02 -70.175  < 2e-16 ***
## DIA_NOMBREJUEVES     4.584e-01  1.007e-02  45.501  < 2e-16 ***
## DIA_NOMBRELUNES      4.244e-01  1.022e-02  41.505  < 2e-16 ***
## DIA_NOMBREMARTES     4.882e-01  1.001e-02  48.779  < 2e-16 ***
## DIA_NOMBREMIERCOLES  4.636e-01  1.005e-02  46.134  < 2e-16 ***
## DIA_NOMBRESABADO     3.979e-01  1.017e-02  39.125  < 2e-16 ***
## DIA_NOMBREVIERNES    5.057e-01  9.973e-03  50.700  < 2e-16 ***
## PERIODO             -6.002e-02  8.733e-03  -6.873 6.28e-12 ***
## ESPECIAL            -2.453e-01  1.020e-02 -24.060  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 199129  on 7246  degrees of freedom
## Residual deviance:  13955  on 7232  degrees of freedom
## AIC: 45706
## 
## Number of Fisher Scoring iterations: 4

Para la evaluacion del modelo se usa nuevamente el MSE y evaluamos para el conjunto de entrenamiento y validación y calculamos la variación porcentual entre los dos:

#MSE para entrenamiento
prediccion_Train <- predict(modelo_poisson, type = "response")


MSE_Train <- mean((prediccion_Train - grupo_1$Accidentes)^2)

# MSE Validacion
grupo_validacion <- Test %>% group_by(FECHA, CLASE, DIA_NOMBRE,
                              PERIODO, ESPECIAL) %>%
  summarise(Accidentes = n()) %>%
  arrange(FECHA)

prediccion_Test <- predict(object = modelo_poisson, newdata = grupo_validacion[,-6], 
                           type = "response")
MSE_Test <- mean((prediccion_Test - grupo_validacion$Accidentes)^2)


print(c("MSE Test" = MSE_Test, "MSE Train" = MSE_Train))
##  MSE Test MSE Train 
##  55.69298  57.50773
Variacion <- ((MSE_Train - MSE_Test)/(MSE_Test*2))*100
print(c("Variacion porcentual" = Variacion))
## Variacion porcentual 
##             1.629249

3.6 Pronosticos

Dado que los modelos son temporales por lo que entre más alejado sea el horizonte de pronóstico menos acertada será la predicción se tomará el año 2019 para predecir con los modelos entrenados con 4 años (2014-2017). Los datos pronosticados se guardarán en un csv que será utilizado para la visualización e interacción del usuario en la aplicación VEHAPP

#====================================
# Prediccion 2019 con modelo poisson
#=====================================
Datos_2019 <- read.csv("Datos/Otros/Predecir_2019.csv", sep = ",", encoding = "UTF-8")
Datos_2019 <- Datos_2019[,-1]
Datos_2019$CLASE <- gsub(pattern = "CAIDA OCUPANTE", replacement = "CAIDA_OCUPANTE", x = Datos_2019$CLASE)

Datos_2019$FECHA <- as.Date(Datos_2019$FECHA, format = "%Y-%m-%d")
Datos_2019$CLASE <- as.factor(Datos_2019$CLASE)
Datos_2019$DIA_NOMBRE <- as.factor(Datos_2019$DIA_NOMBRE)
Datos_2019$PERIODO <- as.integer(Datos_2019$PERIODO)

prediccion2019 <- predict(object = modelo_poisson, newdata = Datos_2019,
                          type = "response")
PREDICCIONES_DIARIAS_2019 <- Datos_2019 %>% 
  mutate(ACCIDENTES = round(prediccion2019,0))


#write.csv(PREDICCIONES_DIARIAS_2019, "Predicciones_diarias_2019.csv")

  1. Conclusiones

En este trabajo se usó información de accidentalidad comprendida de 2014 a 2018, con el fin de buscar un modelo que se ajustara bien a los datos y tuviera la capacidad de predecir años siguientes, en la busqueda se encontró que el model XGboost tiene buena capacidad de predicción en las resoluciones diarias, semanales y mensuales que son agrupadas con conteos y se encontró que el modelo generalizado con distribución Poisson es el mejor modelo para pronosticar dia a dia utilizando fechas.

En la caracterización de los accidentes por barrios se encontró que la diferencia más significativa entre estos es la cantidad de accidentes que ocurren y se llegó a la conclusión que entre más cercano se está a la zona centro más propenso está a sufrir un accidente y que las vias rapidas de la ciudad son centro constante de accidentes también.

  • Los choques son los tipos de accidentes que más se presentan en Medellín
  • El día más grave para la movilidad es el viernes y los días más tranquilos son el domingo y los lunes aunque están influenciados por los festivos
  • Semana Santa se puede considerar como una semana atipica dado que reduce la accidentalidad respecto a las semanas comunes.

Como observación final, el uso de técnicas estadísticas depende mucho de la calidad del dato y de como estos se hayan recolectado y tratado.

  1. Bibliografia

  • James, G., Witten, D., Hastie, T., & Tibshirani, R. (2013). An introduction to statistical learning (Vol. 112, p. 18). New York: springer.
  • Ospina J. (2019) Introducción a la analítica predictiva, Vehículos registrados en el RUNT. Material de clase especialización en Analítica Unal.
  • Concejo de Medellín. La Alta accidentalidad en Medellín se convirtió en un problema de Salud Pública. Boletin 070. 2018
  • IBM Corporation (2015). icesi.edu.co. Obtenido de http://gforge.icesi.edu.co/ASUMDM_External/index.html
  • IBM Corporation. (2016). Analytics Solutions Unified Method. Implementations with agile principles. New York: IBM Analytics. Recuperado en Octubre de 2020, de ftp://ftp.software.ibm.com/software/data/sw-library/services/ASUM.pdf
  • Garcia, H., Vera, Y., Zuluaga, Lina., (2011). Características de los accidentes de tránsito con personas lesionadas atendidas en un hospital de tercer nivel de Medellín, 1999-2008. Bogotá: Rev.Gerenc.polit.Salud. Obtnido de: http://www.scielo.org.co/pdf/rgps/v10n21/v10n21a07.pdf
  • RUNT (2019). RUNT EN CIFRAS. Obtenido de: https://www.runt.com.co/runt-en-cifra

Software R